#include "SpectraSTSearch.hpp"
#include "SpectraSTLibEntry.hpp"
#include "FileUtils.hpp"

#include <algorithm>
#include <math.h>
#include <sstream>


/*

Program       : Spectrast
Author        : Henry Lam <hlam@systemsbiology.org>                                                       
Date          : 03.06.06 


Copyright (C) 2006 Henry Lam

This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.

This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA

Henry Lam
Insitute for Systems Biology
1441 North 34th St. 
Seattle, WA  98103  USA
hlam@systemsbiology.org

*/

/* Class: SpectraSTSearch
 * 
 * Implements the search for one query spectrum.
 *  
 */


using namespace std;


extern bool g_verbose;

// constructor
SpectraSTSearch::SpectraSTSearch(SpectraSTQuery* query,	SpectraSTSearchParams& params, SpectraSTSearchOutput* output) :
  m_query(query),
  m_params(params), 
  m_output(output),
  m_candidates() {
  
  m_candidates.clear();	
  
}

// destructor
SpectraSTSearch::~SpectraSTSearch() {
  if (m_query) delete m_query;
	
  for (vector<SpectraSTCandidate*>::iterator i = m_candidates.begin(); i != m_candidates.end(); i++) {
    delete (*i);	
  }
}

// search - main function to perform one search
void SpectraSTSearch::search(SpectraSTLib* lib) {

  double precursorMz = m_query->getPrecursorMz();

  if (g_verbose) {
    cout << endl;
    cout << "Now searching query: " << m_query->getName() << " (PrecursorMZ = " << precursorMz;
    
    if (m_query->getDefaultCharge() != 0) {
      cout << "; PrecursorCharge = " << m_query->getDefaultCharge();
      for (int ch = 1; ch <= 7; ch++) {
        if (m_query->isPossibleCharge(ch) && ch != m_query->getDefaultCharge()) cout << "," << ch;
      }
    }
    
    cout << ")" << endl;
  }
  cout << "pMz=" << precursorMz  << endl;  
 
  // retrieves all entries from the library within the tolerable m/z range
  vector<SpectraSTLibEntry*> entries;

  if (m_params.searchMode.compare("bonanza") == 0){
    //<DEBUG>
    cout << "Import all entries of seeds!..." ;
    //</DEBUG>
    //lib->getAllEntries(entries);
    double lowMz = precursorMz - m_params.blindModMzTolerance;
    if (lowMz < 0){
      lowMz = 0;
    }
    double highMz = precursorMz + m_params.blindModMzTolerance;
    lib->retrieve(entries, lowMz, highMz);
    //<DEBUG>
    cout << "retrieve candidates in mz range [" << lowMz <<"," << highMz << "]" << endl;
    cout << entries.size() << " seeds are imported!" <<endl;
    //</DEBUG>
  }else{

    //default mode
    double lowMz = precursorMz - m_params.indexRetrievalMzTolerance;
    if (m_params.indexRetrievalUseAverage) lowMz -= 1.0;
    double highMz = precursorMz + m_params.indexRetrievalMzTolerance;
  
    lib->retrieve(entries, lowMz, highMz);
  }  
  if (g_verbose) {
    cout << "\tFound " << entries.size() << " candidate(s)... " << " Comparing... ";
    cout.flush();
  }
  
  // for all retrieved entries, do the necessary filtering, add the good ones to m_candidates
  for (vector<SpectraSTLibEntry*>::iterator i = entries.begin(); i != entries.end(); i++) {
   
    if (!(m_params.ignoreChargeOneLibSpectra && (*i)->getCharge() == 1) &&
	!(!m_params.expectedCysteineMod.empty() && m_params.expectedCysteineMod != "ICAT_cl" && (*i)->isCleavableICAT()) &&
	!(!m_params.expectedCysteineMod.empty() && m_params.expectedCysteineMod != "ICAT_uc" && (*i)->isUncleavableICAT()) &&
	!(!m_params.expectedCysteineMod.empty() && m_params.expectedCysteineMod != "CAM" && (*i)->isCAMCysteine()) &&
	!(m_params.ignoreAbnormalSpectra && (*i)->getStatus() != "Normal") &&
	!(m_params.ignoreSpectraWithUnmodCysteine && (*i)->hasUnmodifiedCysteine()) &&
	(m_params.searchAllCharges || m_query->isPossibleCharge((*i)->getCharge()))) {	
	
      // apply library spectrum simplification
      double retained = (*i)->getPeakList()->simplify(m_params.filterLibMaxPeaksUsed, 999999999);
    }
   
    SpectraSTCandidate* newCandidate = new SpectraSTCandidate(*i, m_params);	
    m_candidates.push_back(newCandidate);
    
  } //loop for entries
  
  
  
  //<DEBUG>
  cout << "Compare query to " << m_candidates.size() << " candidates!" << endl;
  //</DEBUG>
  // compare query to each candidate by calculating the dot product
  for (vector<SpectraSTCandidate*>::iterator i = m_candidates.begin(); i != m_candidates.end(); i++) {

    SpectraSTLibEntry* entry = (*i)->getEntry();
    int charge = entry->getCharge();
    //<DEBUG>
    cout << "libID = " << entry->getLibId() << "\t";
    cout << "fullName = " << entry->getFullName() << "\t";
    //</DEBUG>
    double dot = m_query->getPeakList(charge)->compare(entry->getPeakList(), (*i)->getSimScoresRef(), m_params);
    //<DEBUG>
    cout << "dot=" << dot <<endl;
    //</DEBUG>
    if (m_params.indexRetrievalUseAverage) {
      ((*i)->getSimScoresRef()).precursorMzDiff = precursorMz - entry->getAveragePrecursorMz();          	 
    } else {
      ((*i)->getSimScoresRef()).precursorMzDiff = precursorMz - entry->getPrecursorMz();    
    }
    (*i)->setSortKey(dot);	
    
  }

  // sort the hits by the sort key 
  // (in this case, the value of "dot" returned by the SpectraSTPeakList::compare function)
  sort(m_candidates.begin(), m_candidates.end(), SpectraSTCandidate::sortPtrsDesc);
  
  if (m_params.detectHomologs > 1) {
    detectHomologs();
  } else {
    if (m_candidates.size() > 1) {
      (m_candidates[0]->getSimScoresRef()).firstNonHomolog = 2;
    }
  }
 
  calcHitsStats();
  
  // to save time, we can throw away all except the top N hits
  // this is assuming nothing beyond N will be worth rescuing to the top due to a favorable dot_bias
  // unsigned int numHitsToKeep = m_params.detectHomologs + 6;
  //while (m_candidates.size() > numHitsToKeep) {
  //  delete (m_candidates.back());
  //  m_candidates.pop_back();
  //}
  
  // now that they are sorted, can calculate delta dots.
  calcDeltaSimpleDots();
	
  // sort again by F value (calculated in calcDeltaSimpleDots)
  sort(m_candidates.begin(), m_candidates.end(), SpectraSTCandidate::sortPtrsDesc);
  
  //consider the top one hit
  //check if the dot > 0.4, if less, skip the below steps
  //check if abs(massdiff) > 0.5, if less, skip the below steps
  //add massDiff on each amino acid , generate the modToken and check if any possible match
  //if no matchs found, skip the below steps (or set dot=0 when big massDiff)
  //if one match found, just use it as modification
  //if multiple matches found, for each mod case try to annotate the peaks based this modification and calculate the total intenisity and number of annotated peaks
  //use the highest intensity / most peaks as the top hit,modify the hit name
  if (m_candidates.empty()){
    return;
  }
  
  vector<SpectraSTCandidate*>::iterator itCandidate = m_candidates.begin();
  
  if ((((*itCandidate)->getSimScoresRef()).dot > 0.4)||((*itCandidate)->getEntry()->getCharge()>3)){
    double massDiff = ((*itCandidate)->getSimScoresRef()).precursorMzDiff * (*itCandidate)->getEntry()->getCharge();
    // NOTE: Here, hack code for mztolerance , it should be adjusted when high/low accuracy data used.
    if (fabs(massDiff) > 0.5){
      //<DEBUG>
      cout << "topHit dot = "<<((*itCandidate)->getSimScoresRef()).dot << "\t";
      cout << "massDiff = " << massDiff << endl;
      //</DEBUG>
      //1. get peptide, 
      SpectraSTLibEntry* itEntry = (*itCandidate)->getEntry();
      Peptide* p = itEntry->getPeptidePtr();
      string strip = p->stripped;
      
      //2. get modename
      vector<string> modNames; 
      Peptide::findModNamesByMass(modNames,massDiff,0.5);//new
     
      vector<string> negModNames;
      if (massDiff < 0){
	Peptide::findModNamesByMass(negModNames,fabs(massDiff),0.5);
      }
      
      map<unsigned int, vector<string> > modTokenCandidates;
      map<unsigned int, vector<string> >::iterator itModTokenCandidates;
      unsigned int modTokenCandidatesNum = 0;
      bool isUnknownModType = false;

      if ((!modNames.empty())||(!negModNames.empty())){
	//known modification can be found!
	//3. foreach amino acid, get original mass ,current mass (check if it is already modified), calculate the new modified mass , generate new mod token and see if it can be found in the map of all modfication token list
	//4. switch: nothing matched, one matched, multiply matched
	for (unsigned int i = 0; i < strip.length() ; i++){
	  //find a existed modified position, impossible to add extra modification in the same place
	  if ( (p->isModsSet)&& (p->mods.find(int(i)) != p->mods.end()) && (!negModNames.empty()) ){
	    // if the amino acid is already modified, just check if it is demodified
	    if (massDiff < 0){
	      string modType =  p->mods.find(int(i))->second;
	      for (vector<string>::iterator itNeg = negModNames.begin(); itNeg != negModNames.end(); itNeg ++){
		
		if (itNeg->compare(modType) == 0){
		  string modToken = strip.substr(i,1);
		  vector<string> modTokenCandidatesAti (1,modToken);
		  modTokenCandidates.insert(pair<unsigned int, vector<string> >( i,modTokenCandidatesAti));
		  modTokenCandidatesNum ++;
		}
	      }
	      
	    }
	    continue;
	  }

	  // if no reasonable mod names can be found, try to search next amino acd
	  if (modNames.empty()){
	    continue;
	  }

	  char aa = strip[i];
	  double aaMass = 0;
	  map<char,double>::iterator aaFound = Peptide::AAMonoisotopicMassTable->find(aa);
	  if (aaFound != Peptide::AAMonoisotopicMassTable->end()){
	    //get mass of current amino acid
	    aaMass = aaFound->second;
	    //look for known reasonable modToken candidates
	    for(vector<string>::iterator itModName = modNames.begin(); itModName != modNames.end(); itModName++){
	      string modToken = p->getModToken(aa,(*itModName));
	      //cout << "modName=" << (*itModName) << endl;
	      //cout << "modToken=" << modToken << endl;
	      map<string,string>::iterator modTokenFound = Peptide::modTokenTable->find(modToken);
	      if (modTokenFound != Peptide::modTokenTable->end()){
		
		itModTokenCandidates = modTokenCandidates.find(i);
		if (itModTokenCandidates != modTokenCandidates.end() ){
		  //cout << "same site add modtoken=" << modToken << endl;
		  // in the same amino acid site, we already have a mod candidate, we need to add a new candidate
		  itModTokenCandidates->second.push_back(modToken);
		}else{
		  //cout << "new site add modtoken=" << modToken << endl;
		  // add a new candidate in this site
		  vector<string> modTokenCandidatesAti;
		  modTokenCandidatesAti.push_back(modToken);
		  modTokenCandidates.insert(pair<unsigned int,vector<string> >(i, modTokenCandidatesAti));
		}
		//the new modtoken has been added to modTokenCandidates
		modTokenCandidatesNum ++;
	      }// if
	    }// for loop of modnames
	  }else{
	    cout << "ERROR: unknown amino acid in the peptide!" << endl;
	  }// if amino acid is legal
	}// for loop of stripped sequence
      }// if modnames are not empty
      
      if ( modTokenCandidatesNum == 0 ){
	isUnknownModType = true;
	// modnames are empty, we can not find any known modification
	// this may be a new kind of modification, try to find a suitable site for it
	// add modToken for each amino acid ???
	for (unsigned int i = 0; i < strip.length() ; i++){
	  //find a modified position, impossible to add modification here
	  if ((p->isModsSet)&& (p->mods.find(int(i)) != p->mods.end())){
	    continue;
	  }
	  char aa = strip[i];
	  double aaMass = 0;
	  map<char,double>::iterator aaFound = Peptide::AAMonoisotopicMassTable->find(aa);
	  if (aaFound != Peptide::AAMonoisotopicMassTable->end()){	  
	    //get mass of current amino acid
	    aaMass = aaFound->second;
	    //give a legal modToken format
	    stringstream ssMassDiff;
	    ssMassDiff << massDiff;
	    string modToken = p->getModToken(aa,ssMassDiff.str(),true);
	    vector<string> modTokenCandidatesAti;
	    modTokenCandidatesAti.push_back(modToken);
	    modTokenCandidates.insert(pair<unsigned int,vector<string> >(i, modTokenCandidatesAti));	    
	    modTokenCandidatesNum ++;
	    stringstream ssModType;
	    ssModType << i << "_plus_" << massDiff;
	    Peptide::addModTokenToTables(modToken,ssModType.str() ,true, massDiff);
	    Peptide::modTablesChanged = true;
	  }	
	}
      }
	
      //compare all the modification candidates and see which one is best
      if (!modTokenCandidates.empty()){
	//<DEBUG>
	cout << "modToken Candidates List: " << endl;
	for(itModTokenCandidates = modTokenCandidates.begin(); itModTokenCandidates != modTokenCandidates.end(); itModTokenCandidates ++){
	  unsigned int pos = itModTokenCandidates->first;
	  vector<string>::iterator itModTokensAti;
	  for (itModTokensAti = itModTokenCandidates->second.begin(); itModTokensAti != itModTokenCandidates->second.end(); itModTokensAti++){
	    cout << "pos=" << pos << "\t" << (*itModTokensAti) << endl;
	  }
	}
	cout << "END: modToken Candidates List " << endl;
	//</DEBUG>
	vector< pair<string,double> > scores;//pair<"pos,modToken",score> 
	//loop in all amino acids
	for(itModTokenCandidates = modTokenCandidates.begin(); itModTokenCandidates != modTokenCandidates.end(); itModTokenCandidates ++ ){
	  //loop for all mod candidates in one site
	  vector<string>::iterator i;
	  for (i = itModTokenCandidates->second.begin(); i !=itModTokenCandidates->second.end(); i++){
	    Peptide* newPep = new  Peptide( *( itEntry->getPeptidePtr() ) );
	    unsigned int pos = itModTokenCandidates->first;
	    if (!newPep->addMod(pos, (*i))){
	      continue;
	    }
	    m_query->getPeakList()->setPeptidePtr(newPep);
	    m_query->getPeakList()->annotate(true,false);
	    double annotationScore =   m_query->getPeakList()->getAnnotationScore();
	    stringstream ssMod;
	    ssMod << pos << "," << (*i); 
	    scores.push_back(pair<string,double>(ssMod.str(),annotationScore));
	    delete newPep;
	  }
	}
	
	if (!scores.empty()){ 
	  // choose the most suitable annotation
	  sort(scores.begin(),scores.end(), SpectraSTPeakList::sortByAnnotationScoreDesc);
	  //<DEBUG>
	  cout << "The best modToken candidate is" << (scores.begin())->first << "\t" << "The annoatation score is  " << (scores.begin())->second << endl;
	  //</DEBUG>
	  
	  // assign the correct peptide to simScore and update precursorMzDiff
	  // annotation score must > threshold for an unknown modification, now just use zero, hack code
	  vector< pair<string,double> >::iterator itScore = scores.begin();
	  if ((itScore->second < 0)&&(isUnknownModType)){
	    //consider to redo a regular search
	    cout << "WARN: too low annotation score=" << itScore->second << endl;
	 
	  }else{
	    size_t site = 0;
	    unsigned int pos = atoi(nextToken(itScore->first,site,site,",").c_str());
	    string modToken = nextToken(itScore->first,site+1,site);
	    Peptide* newPep = new Peptide (*( itEntry->getPeptidePtr() ));
	    newPep->addMod(pos,modToken);
	    (*itCandidate)->getSimScoresRef().newPep = newPep;
	    (*itCandidate)->getSimScoresRef().precursorMzDiff = m_query->getPrecursorMz() - newPep->monoisotopicMZ();
	    //<DEBUG>
	    //cout << "query_precursorMz=" << m_query->getPrecursorMz() << "\t";
	    //cout << "newPep_mz=" << newPep->monoisotopicMZ() << endl;
	    //</DEBUG>
	  }	  
	}else{
	  //consider to redo a regular search as it can not be a modified version of the current libentry
	  cout << "WARN: no modfications with annotatoin scores found!" << endl;
	}
      }else{
	//consider to redo a regular search by using bonanza score without partial-match in a regular mz tolerance
	cout << "WARN: no modtoken candidates found!" << endl;
      }
    }//if massDiff > 0.5
  } else{
    //<DEBUG>
    cout << "dot =" << ((*itCandidate)->getSimScoresRef()).dot << endl;
    //</DEBUG>
  }


}

// calcDeltaSimpleDots - calculates the delta dots
void SpectraSTSearch::calcDeltaSimpleDots() {

  if (m_candidates.empty()) {
    return;
  }
  
  for (unsigned int rank = 0; rank < (unsigned int)(m_candidates.size()); rank++) {
    SpectraSTSimScores& s = m_candidates[rank]->getSimScoresRef();
    
    if (s.firstNonHomolog > 0) {
      // delta is dot - dot(first non-homologous hit lower than this one)
      s.delta = s.dot - (m_candidates[s.firstNonHomolog - 1]->getSimScoresRef()).dot;
    } else {
      s.delta = 0.0;
    }
    
    double fval = s.calcFval(m_params.fvalFractionDelta);
    
    m_candidates[rank]->setSortKey(fval);
  }
}

// calcHitsStats - calculates such things as the mean and stdev of the dots of all candidates
void SpectraSTSearch::calcHitsStats() {
  
  unsigned int numHits;
  double mean;
  double stdev;
  
  if (m_candidates.empty()) {
    numHits = 0;
    mean = 0.0;
    stdev = 0.0;
  } else {
    
    double totalDot = 0;
    double totalSqDot = 0;
    for (vector<SpectraSTCandidate*>::iterator i = m_candidates.begin(); i != m_candidates.end(); i++) {
      SpectraSTSimScores& s = (*i)->getSimScoresRef();
      totalDot += s.dot;
      totalSqDot += s.dot * s.dot;
    }
    numHits = (unsigned int)(m_candidates.size());
    mean = totalDot / (double)numHits;
    stdev = totalSqDot / (double)numHits - mean * mean;
    if (stdev > 0.000001) stdev = sqrt(stdev);
  }
  for (vector<SpectraSTCandidate*>::iterator i = m_candidates.begin(); i != m_candidates.end(); i++) {
    SpectraSTSimScores& s = (*i)->getSimScoresRef();
    s.hitsNum = numHits;
    s.hitsMean = mean;
    s.hitsStDev = stdev;
    
  }
  	
}

// detectHomologs - try to see if the lower hits are homologous (or identical) to the first one
void SpectraSTSearch::detectHomologs() {
  
  if (m_candidates.empty() || m_candidates.size() == 1) return;
	
  Peptide* topHit = m_candidates[0]->getEntry()->getPeptidePtr();
  
  if (!topHit) {
    // not a peptide. no notion of homology
    (m_candidates[0]->getSimScoresRef()).firstNonHomolog = 2;  
    return;
  }
    
  bool homologFound = false;
  unsigned int curRank = 0;
  
  // go down the hit lists until hitting a nonhomologous hit
  do {
    homologFound = false;
    curRank++;
    
    Peptide* thisHit = m_candidates[curRank]->getEntry()->getPeptidePtr();
    
    if (!thisHit) {
      // not a peptide. definitely nonhomologous
      break;
    }
    
    int identity = 0;
    if (*topHit == *thisHit) {
      // identical!
      homologFound = true;
    } else if ((topHit->stripped.length() > thisHit->stripped.length() && topHit->isSubsequence(*thisHit, true)) ||
	       (thisHit->stripped.length() < topHit->stripped.length() && thisHit->isSubsequence(*topHit, true))) {
      // one is subsequence of the other!
      homologFound = true;
    } else if (topHit->isHomolog(*thisHit, 0.7, identity)) {
      homologFound = true;
      
    } 
    
  } while (homologFound && curRank < m_params.detectHomologs - 1 && curRank < (unsigned int)(m_candidates.size()) - 1);
  
  // setting the field firstNonHomolog for all the homologs found
  for (unsigned int rank = 0; rank < curRank; rank++) {
    (m_candidates[rank]->getSimScoresRef()).firstNonHomolog = curRank + 1;
  }
			
}



// print - prints out the search result
void SpectraSTSearch::print() {
  
	
  if (m_candidates.empty() || (!(m_candidates[0]->passTopHitFvalThreshold()))) {
    // no hit (or no hit above threshold
    if (g_verbose) {
      cout << " DONE! Top hit: NO_MATCH";
      cout.flush();
    }
    if (m_params.hitListExcludeNoMatch) {
      // told to exclude all NO_MATCH's from the final output, so just return
      return;
    }
  }
  
  // set the assumedCharge to that of the top hit, if any
  int assumedCharge = 0;
  if (!m_candidates.empty() && m_candidates[0]->passTopHitFvalThreshold()) {
    assumedCharge = m_candidates[0]->getEntry()->getCharge();
  }
  
  
  string name = m_query->getName();
  double precursorMz = m_query->getPrecursorMz();
  double rt = m_query->getRetentionTime();
  // prints all the query information
  m_output->printStartQuery(name, precursorMz, assumedCharge, rt);
  
  
  if (m_candidates.empty() || !(m_candidates[0]->passTopHitFvalThreshold())) {
    // no hit or no hit above threshold
    SpectraSTSimScores emptyScore;
		
    m_output->printHit(name, 1, NULL, emptyScore);
		
  } else {
    // print the top hit!
    if (g_verbose) {
      cout << " DONE! Top hit: " << m_candidates[0]->getEntry()->getFullName() << " (Dot = " << m_candidates[0]->getSortKey() << ")";
      cout.flush();
    }
    
    m_output->printHit(name, 1, m_candidates[0]->getEntry(), m_candidates[0]->getSimScoresRef());
    
    /*    
    if (m_params.saveSpectra) {
      // also save the query spectrum and the spectrum for the top hit
      string querySpectrumFileName(m_output->getOutputPath() + name + ".0.msp");
      m_query->printQuerySpectrum(querySpectrumFileName);
      
      string topMatchSpectrumFileName(m_output->getOutputPath() + name + ".1.msp");		
      m_candidates[0]->printSpectrum(topMatchSpectrumFileName);
    }
    */    
    
    if (!m_params.hitListOnlyTopHit) { 
      // told to print the lower hits too
      
      unsigned int firstNonHomolog = (m_candidates[0]->getSimScoresRef()).firstNonHomolog;
      
      for (unsigned int rank = 2; rank <= (unsigned int)(m_candidates.size()); rank++) {		
        if (m_candidates[rank - 1]->passLowerHitsFvalThreshold() || (m_params.hitListShowHomologs && rank < firstNonHomolog)) {
          m_output->printHit(name, rank, m_candidates[rank - 1]->getEntry(), m_candidates[rank - 1]->getSimScoresRef());
        } else {
          // below threshold now, stop printing
          break;
        }				
      }
    }
    
  }
  
  // print the closing tags to finish up
  m_output->printEndQuery(name);
  
  
}

// isLikelyGood. Simply returns if the F value is above 0.5, indicating a likely good hit (not always!)
bool SpectraSTSearch::isLikelyGood() {
  if (m_candidates.size() > 0 && m_candidates[0]->getSortKey() >= 0.5) {
    return (true);
  } else {
    return (false);
  }
  return (false);
}

// isLikelyBad. Simply returns if the F value is below 0.2.
bool SpectraSTSearch::isLikelyBad() {
  if (m_candidates.size() > 0 && m_candidates[0]->getSortKey() < 0.2) {
    return (true);
  } else {
    return (false);
  }
  return (false);
}

// isDecoy. returns if the hit is a decoy entry
bool SpectraSTSearch::isDecoy(unsigned int rank) {
  if ((unsigned int)(m_candidates.size()) >= rank) {
    string spec("");
    string rem("");
    return ((m_candidates[rank - 1]->getEntry()->getOneComment("Spec", spec) && spec == "Decoy") ||
	    (m_candidates[rank - 1]->getEntry()->getOneComment("Remark", rem) && rem.substr(0,5) == "DECOY"));
  } else {
    return (false);
  }
  return (false);
}

bool SpectraSTSearch::isSingleton(unsigned int rank) {
  if ((unsigned int)(m_candidates.size()) >= rank) {
    string nreps("");
    return (m_candidates[rank - 1]->getEntry()->getNrepsUsed() == 1);
  } 
  return (false);
}
